In [1]:
import cvxpy as cvx
import numpy as np
import matplotlib.pylab as plt
from typing import Tuple
import altair as alt
from vega_datasets import data
import pandas as pd
import plotly.graph_objects as go

Maximum likelihood estimation of an increasing nonnegative signal

In [2]:
# create problem data
N = 100
# create an increasing input signal
xtrue = np.zeros(N)
xtrue[0:40] = 0.1
xtrue[50] = 2
xtrue[69:80] = 0.15
xtrue[80] = 1
xtrue = np.cumsum(xtrue)

# pass the increasing input through a moving-average filter
# and add Gaussian noise
h = np.array([1, -0.85, 0.7, -0.3])
k = len(h)
yhat = np.convolve(h, xtrue)
In [3]:
# Data ported from matlab file
y = yhat[0:-3] + np.array(
    [
        -0.43,
        -1.7,
        0.13,
        0.29,
        -1.1,
        1.2,
        1.2,
        -0.038,
        0.33,
        0.17,
        -0.19,
        0.73,
        -0.59,
        2.2,
        -0.14,
        0.11,
        1.1,
        0.059,
        -0.096,
        -0.83,
        0.29,
        -1.3,
        0.71,
        1.6,
        -0.69,
        0.86,
        1.3,
        -1.6,
        -1.4,
        0.57,
        -0.4,
        0.69,
        0.82,
        0.71,
        1.3,
        0.67,
        1.2,
        -1.2,
        -0.02,
        -0.16,
        -1.6,
        0.26,
        -1.1,
        1.4,
        -0.81,
        0.53,
        0.22,
        -0.92,
        -2.2,
        -0.059,
        -1,
        0.61,
        0.51,
        1.7,
        0.59,
        -0.64,
        0.38,
        -1,
        -0.02,
        -0.048,
        4.3e-05,
        -0.32,
        1.1,
        -1.9,
        0.43,
        0.9,
        0.73,
        0.58,
        0.04,
        0.68,
        0.57,
        -0.26,
        -0.38,
        -0.3,
        -1.5,
        -0.23,
        0.12,
        0.31,
        1.4,
        -0.35,
        0.62,
        0.8,
        0.94,
        -0.99,
        0.21,
        0.24,
        -1,
        -0.74,
        1.1,
        -0.13,
        0.39,
        0.088,
        -0.64,
        -0.56,
        0.44,
        -0.95,
        0.78,
        0.57,
        -0.82,
        -0.27,
    ]
)
In [4]:
xtrue
Out[4]:
array([0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.  , 1.1 ,
       1.2 , 1.3 , 1.4 , 1.5 , 1.6 , 1.7 , 1.8 , 1.9 , 2.  , 2.1 , 2.2 ,
       2.3 , 2.4 , 2.5 , 2.6 , 2.7 , 2.8 , 2.9 , 3.  , 3.1 , 3.2 , 3.3 ,
       3.4 , 3.5 , 3.6 , 3.7 , 3.8 , 3.9 , 4.  , 4.  , 4.  , 4.  , 4.  ,
       4.  , 4.  , 4.  , 4.  , 4.  , 4.  , 6.  , 6.  , 6.  , 6.  , 6.  ,
       6.  , 6.  , 6.  , 6.  , 6.  , 6.  , 6.  , 6.  , 6.  , 6.  , 6.  ,
       6.  , 6.  , 6.  , 6.15, 6.3 , 6.45, 6.6 , 6.75, 6.9 , 7.05, 7.2 ,
       7.35, 7.5 , 7.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65,
       8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65, 8.65,
       8.65])

Formulate the problem of finding the maximum likelihood estimate of x, given y, taking into account the prior assumption that x is nonnegative and monotonically nondecreasing, as a convex optimization problem.

$ \begin{align} \mathbf{maximize} \; &\, \; p(x) = -(m/2) \cdot log(2 \pi \sigma^2) -\frac{1}{2\sigma^2} \biggl\lvert\biggl\lvert \; \sum_{\tau=1}^k {h(\tau) x(t-\tau)} -y \;\biggr\rvert\biggr\rvert_2 \newline \mathbf{subject \; to}: \, & \; x_i>=0 \newline & x_{i+1} > x_i \end{align} $

This is equivalent to:

$ \begin{align} \mathbf{minimize} \; &\, \; \hat{y} -y \newline \mathbf{subject \; to}: \, & \; x_i>=0 \newline & x_{i+1} > x_i \end{align} $

where $\hat{y} = \biggl\lvert\biggl\lvert \; \sum_{\tau=1}^k {h(\tau) x(t-\tau)}\;\biggr\rvert\biggr\rvert_2 $

In [5]:
# Define variables
x_constrained = cvx.Variable(N)
x_free = cvx.Variable(N)
In [6]:
# Define objectives
obj_constrained = cvx.Minimize(
    cvx.norm(cvx.conv(h, x_constrained)[0:-3] - y[:, np.newaxis])
)
obj_free = cvx.Minimize(cvx.norm(cvx.conv(h, x_free)[0:-3] - y[:, np.newaxis]))
In [7]:
# Define constraints
constraints = [
    x_constrained >= 0,
    x_constrained[1:] >= x_constrained[:-1]
]
In [8]:
# Form and solve the constrained problem
prob_constrained = cvx.Problem(obj_constrained, constraints)
prob_constrained.solve()  # Returns the optimal value.
Out[8]:
7.627761550201712
In [9]:
# Form and solve the free problem
prob_free = cvx.Problem(obj_free)
prob_free.solve()  # Returns the optimal value.
Out[9]:
-4.8647294852831714e-12
In [10]:
print("optimal x_constrained", x_constrained.value)
print("optimal x_free", x_free.value)
optimal x_constrained [3.23379604e-09 3.92315295e-09 8.12553256e-09 6.63452349e-01
 6.63452350e-01 8.35625772e-01 1.22371276e+00 1.22371277e+00
 1.22371277e+00 1.22371277e+00 1.22371277e+00 1.66339981e+00
 1.66339981e+00 2.18310707e+00 2.18310707e+00 2.18310707e+00
 2.18310707e+00 2.18310707e+00 2.18310707e+00 2.18310707e+00
 2.18310709e+00 2.18310709e+00 2.18310709e+00 2.94762867e+00
 2.94762867e+00 2.94762867e+00 2.94762867e+00 2.94762867e+00
 2.94762867e+00 3.50000263e+00 3.50000264e+00 3.52677431e+00
 3.71811884e+00 3.71811884e+00 3.71811884e+00 3.71811884e+00
 3.71811885e+00 3.71811885e+00 3.71811885e+00 3.71811885e+00
 3.71811885e+00 3.71811885e+00 3.71811885e+00 3.90951969e+00
 3.90951969e+00 3.90951969e+00 3.90951969e+00 3.90951970e+00
 3.90951970e+00 3.90951971e+00 5.21393611e+00 5.95206795e+00
 5.95206799e+00 5.95206800e+00 5.95206800e+00 5.95206800e+00
 5.95206800e+00 5.95206800e+00 5.95206800e+00 5.95206800e+00
 5.95206800e+00 5.95206800e+00 6.56898271e+00 6.56898271e+00
 6.56898271e+00 6.60114286e+00 6.60114286e+00 6.60114286e+00
 6.60114287e+00 6.60114287e+00 6.60114287e+00 6.60114287e+00
 6.60114287e+00 6.60114287e+00 6.60114287e+00 6.78681293e+00
 7.04078786e+00 7.06911523e+00 8.42123213e+00 8.42123213e+00
 8.57656152e+00 8.57656152e+00 8.57656154e+00 8.57656154e+00
 8.57656154e+00 8.57656154e+00 8.57656154e+00 8.57656154e+00
 8.68563140e+00 8.68563140e+00 8.68563140e+00 8.68563140e+00
 8.68563140e+00 8.68563141e+00 8.91295278e+00 8.91295278e+00
 8.91295279e+00 8.91295280e+00 8.91295280e+00 8.91295280e+00]
optimal x_free [-0.33       -1.8655     -1.024675    0.88087625  0.11636731  0.73989634
  2.43171765  2.02094275  1.09756789  1.00278807  1.14035516  2.02162061
  1.38096532  3.10579265  2.99973421  1.81500882  2.44468135  2.79139323
  2.18990995  0.9458526   1.58845571  1.29006352  2.27839077  4.46512441
  3.30750126  2.69430622  4.13444673  2.59551574  0.35036754  1.97628541
  3.55324004  4.22696451  4.36853743  4.43535367  5.24016377  5.34495287
  5.4757014   3.40992831  2.82893395  3.73535446  2.80777602  3.11054168
  2.89912355  4.72920865  4.31360337  3.95585388  4.18171603  3.35944192
  1.11508058  1.99172395  5.12024153  6.40252271  6.86549242  7.88997512
  7.71139097  5.75134747  5.5376642   5.29448864  5.62935464  5.99210866
  6.0411337   5.60929398  6.73673889  5.01206238  4.957324    6.9263034
  7.9508498   7.27700715  6.03775224  6.55343933  7.56959901  6.9980773
  5.91817819  5.86767706  5.29172395  5.86604488  7.1747345   7.98481009
  9.1020879   8.20982799  8.8148353   9.67885677 10.51709195  8.576279
  7.79902982  8.77840763  8.33270931  7.29512652  8.85898331  9.55086005
  9.25298068  8.68262654  7.8859041   7.59857411  8.74094309  8.28407098
  8.7398724   9.57982478  8.44766168  7.5840968 ]
In [11]:
# Plot results
df_true = pd.DataFrame({"t": range(100), "x": xtrue, "type": "true"})
df_constrained = pd.DataFrame(
    {"t": range(100), "x": x_constrained.value, "type": "constrained"}
)
df_free = pd.DataFrame({"t": range(100), "x": x_free.value, "type": "free"})
df = pd.concat([df_true, df_constrained, df_free])


alt.Chart(df).mark_line().encode(
    x="t:T",
    y="x:Q",
    color=alt.Color(
        "type",
        scale=alt.Scale(
            domain=["true", "constrained", "free"], range=["red", "green", "blue"]
        ),
    ),
).properties(width=800, height=450)
Out[11]:

Worst-case probability of loss

Two investments are made, with random returns $R_1$ and $R_2$. Find maximum $p_\mathrm{loss} = prob(R_1 + R_2 <= 0)$.

$R_1$ and $R_2$ have Gaussian marginal distributions with known means $\mu_1$ and $\mu_2$ and known standard deviations $\sigma_1$ and $\sigma_2$.

$R_1$ and $R_2$ are correlated with correllation coefficient $\rho$:

\begin{align} E(R_1 - \mu_1)(R_2 - \mu_2) = \rho \sigma_1 \sigma_2 \end{align}

Check joint Gaussian distribution first:

In [12]:
mu_1 = 8
mu_2 = 20
sigma_1 = 6
sigma_2 = 17.5
rho = -0.25
In [13]:
from scipy.stats import norm
In [14]:
p_loss = norm.cdf(
    0,
    mu_1 + mu_2,
    np.sqrt(sigma_1 * sigma_1 + sigma_2 * sigma_2 + 2 * rho * sigma_1 * sigma_2),
)
p_loss
Out[14]:
0.0499925581431758

Now: Do not assume a joint gaussian distribution. The joint distribution can be anything as long as the marginal distributions are Gaussian:

$$ \begin{align} \mathbf{maximize} \; &\, \; \sum_{R_1 + R_2 <= 0} P_{ij} \\ \mathbf{subject \; to}: \, & \; P_{ij}>=0 \;\;\;\;\;\; \mathrm{(this\ is\ a\ pdf)} \\ & \mathbf(P) * \mathbf{1} = p^{(1)} \;\;\;\;\;\; \mathrm{(marginal\ distribution\ for\ p^{(1)}\ is\ Gaussian)} \\ & \mathbf(P)^T * \mathbf{1} = p^{(2)} \;\;\;\;\;\; \mathrm{(marginal\ distribution\ for\ p^{(2)}\ is\ Gaussian)} \\ & (r - \mu_1 \mathbf{1})^T P (r - \mu_2\mathbf{1}) = \rho \sigma_1 \sigma_2 \;\;\;\;\;\; \mathrm{(correlation)} \end{align} $$

Assuming a discretization of $R_1$ and $R_2$ between $r_1 = -30$ and $r_n = 70$, the marginal distributions for $p^{(1)}$ and $p^{(2)}$ are given by:

$$ \begin{align} p_i^{(1)} = \mathrm{prob}(R_1 = r_i) &=& \frac{\exp( -(r_i-\mu_1)^2 / (2\sigma_1^2))} {\sum_{j=1}^n\exp(-(r_j-\mu_1)^2/(2\sigma_1^2))}\\\,\\ p_i^{(1)} = \mathrm{prob}(R_1 = r_i) &=& \frac{\exp( -(r_i-\mu_2)^2 / (2\sigma_2^2))} {\sum_{j=1}^n\exp(-(r_j-\mu_2)^2/(2\sigma_2^2))} \end{align} $$
In [15]:
# Compute linear grid for approximate solution
m = 100
r = np.linspace(-30, 70, m)
In [16]:
# Marginal probabilities:
p_1 = np.exp(-(r - mu_1) * (r - mu_1) / 2 / sigma_1 / sigma_1) / np.sum(
    np.exp(-(r - mu_1) * (r - mu_1) / 2 / sigma_1 / sigma_1)
)
p_2 = np.exp(-(r - mu_2) * (r - mu_2) / 2 / sigma_2 / sigma_2) / np.sum(
    np.exp(-(r - mu_2) * (r - mu_2) / 2 / sigma_2 / sigma_2)
)
In [17]:
# Define loss mask as indicator matrix for entries that denote loss situations (i.e. r1 + r2 <= 0)
loss_mask = np.sum(np.meshgrid(r, r), axis=0) <= 0
In [18]:
# Define cvx problem
P = cvx.Variable((m, m))
objective = cvx.Maximize(cvx.sum(cvx.multiply(P, loss_mask)))
constraints = [
    P >= 0,
    P @ np.ones(m) == p_1,
    P.T @ np.ones(m) == p_2,
    (r - mu_1) @ P @ (r - mu_2) == rho * sigma_1 * sigma_2,
]

problem = cvx.Problem(objective, constraints)
In [19]:
# Solve problem
problem.solve()
Out[19]:
0.19203644871853245
In [20]:
# Define points for plots
x = np.repeat(r, m)
y = np.tile(r, m)
z = P.value.flatten()
In [21]:
# Create 3D Mesh plot
fig = go.Figure(data=[go.Mesh3d(x=x, y=y, z=z, color="red", opacity=0.50)])
fig.show()
In [22]:
# Create contour plot
fig = go.Figure(data=go.Contour(x=x, y=y, z=z))
fig.update_layout(
    autosize=False,
    width=800,
    height=800,
)
fig.show()
In [ ]: